load("defenseJan2017.Rdata")
defense<-replicate.def

#######
###Installing Necessary Packages
#####

install.packages(c("logistf", "brglm", "mvtnorm","xtable"))

###################
###MODEL PAPER-TABLE 1
#############

library(logistf)
summary(model1<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+
 oecd+labor_full+wcombat,data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))
 
 #############
 ### Generating Results Table
 ##############
 
names(model1)
list.bb<-model1[c(3,1,4,17)]
mat.bb<-matrix(NA,length(list.bb$terms),4)
mat.bb[,1]<-round(list.bb$coefficients,4)
mat.bb[,2]<-round(sqrt(diag(list.bb$var)),4)
mat.bb[,3]<-round(list.bb$prob,4)
mat.bb[,4]<-round(exp(list.bb$coefficients),4)
rownames(mat.bb)<-list.bb$terms
colnames(mat.bb)<-c("Estimate","Std. Error","p-value","exp(estimate)")
library(xtable)
xtable(mat.bb)
 
 
 ########
 ### Relative odds
 ############

#Fatal dispute
exp(-1.94654235)

#Mil Dict.
exp(-1.49101785)

#Female MPs (10% Change in Seatshare)
exp(10* 0.03840136)

#Female Chief Exec
exp(1.11639852)

#################
###PREDICTED VALUES
###########

library(brglm)
model1.br<-brglm(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+
 oecd+labor_full+wcombat,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE)


###
##FUNCTIONS FOR INTERPRETATION WITH PREDICTED VALUES 
###

#loading mvtnorm for plotting
library(mvtnorm)

#computes empirical hpd
emp.hpd<-function(x, conf){
    conf <- min(conf, 1 - conf)
    n <- length(x)
    nn <- round(n * conf)
    x <- sort(x)
    xx <- x[(n - nn + 1):n] - x[1:nn]
    m <- min(xx)
    nnn <- which(xx == m)[1]
    return(c(x[nnn], x[n - nn + nnn]))
}

pred.boot.hazard<-function(model,x.mat,x.values,conf.level,unlink,num.draws,add=FALSE,to.unlink=TRUE,ylim,col,...){
    #model is an lm object
    #x.mat is a matrix of covariate values (gotten from setx for example)
    #where one value changes and all others are held constant
    #x.values is the vector of values for the covariate that is changing
    #conf.level is confidence level
    #unlink is a function to unlink the link function
    #num.draws is the number of random draws
    if(missing(conf.level)){conf.level<-.95}
    if(missing(unlink)){unlink<-function(x){exp(x)}}
    if(to.unlink==FALSE){unlink<-function(x){x}}
    if(missing(num.draws)){num.draws<-10000}
   	beta.hat<-coef(model)
	beta.cov<-vcov(model)
    beta<-rmvnorm(num.draws,beta.hat,beta.cov)
    pred<-unlink(beta%*%t(x.mat))
    hpd.pred<-apply(pred,2,emp.hpd,conf=conf.level)
    med<-apply(pred,2,median)
    ymin<-min(cbind(med,t(hpd.pred)))
    ymax<-max(cbind(med,t(hpd.pred)))
    if(missing(ylim)){ylim=c(ymin,1)}
    if(missing(col)){col="black"}
    if(add==FALSE){
	    plot(x.values,med,ylim=ylim,col=col,...)
	}else{
		points(x.values,med,ylim=ylim,col=col,...)
	}
    points(x.values,hpd.pred[1,],type="l",lty=2,col=col)
    points(x.values,hpd.pred[2,],type="l",lty=2,col=col)
    return(list(x.values=x.values,median.pred=med,hpd.pred=hpd.pred))
}

diff.boot.hazard<-function(model,x.1,x.2,conf.level,unlink,num.draws){
    #model is an lm object
    #x.1 is one vector of covariate values (gotten from setx for example)
    #x.2 is one vector of covariate values (gotten from setx for example)
    #conf.level is confidence level
    #unlink is a function to unlink the link function
    #num.draws is the number of random draws
    if(missing(conf.level)){conf.level<-.95}
    if(missing(unlink)){unlink<-function(x){1-exp(-exp(x))}}
    if(missing(num.draws)){num.draws<-100000}
   	beta.hat<-coef(model)
   	ind.na<-is.na(beta.hat)
   	beta.hat<-beta.hat[!ind.na]
   	x.1<-x.1[!ind.na]
   	x.2<-x.2[!ind.na]
    beta.cov<-vcov(model)
    beta<-rmvnorm(num.draws,beta.hat,beta.cov)
    pred.1<-unlink(beta%*%x.1)
    pred.2<-unlink(beta%*%x.2)
    median.1<-median(pred.1)
    median.2<-median(pred.2)
    d<-pred.2-pred.1 #difference of predicted values
    hpd.diff<-emp.hpd(d,conf.level)
    med.diff<-median(d)
    return(list(median.1=median.1,median.2=median.2,median.diff=med.diff,hpd.diff=hpd.diff))
}

db.2<-function(x.1,model,x.2,conf.level,unlink,num.draws){
    #model is an lm object
    #x.1 is one vector of covariate values (gotten from setx for example)
    #x.2 is one vector of covariate values (gotten from setx for example)
    #conf.level is confidence level
    #unlink is a function to unlink the link function
    #num.draws is the number of random draws
    if(missing(conf.level)){conf.level<-.9}
    if(missing(unlink)){unlink<-function(x){1/(1+exp(-x))}}#Logit Prob
    #if(missing(unlink)){unlink<-function(x){exp(x)}}#Logit Odds
       # if(missing(unlink)){unlink<-function(x){exp(x)}}#C Log Log
    if(missing(num.draws)){num.draws<-100000}
   	beta.hat<-coef(model)
   	ind.na<-is.na(beta.hat)
   	beta.hat<-beta.hat[!ind.na]
   	x.1<-x.1[!ind.na]
   	x.2<-x.2[!ind.na]
    beta.cov<-vcov(model)
    beta<-rmvnorm(num.draws,beta.hat,beta.cov)
    pred.1<-unlink(beta%*%x.1)
    pred.2<-unlink(beta%*%x.2)
    median.1<-median(pred.1)
    median.2<-median(pred.2)
    d<-pred.1/pred.2 #cox hazard ratio
    hpd.diff<-emp.hpd(d,conf.level)
    med.diff<-median(d)
quantile.diff<-quantile(d,c((1-conf.level)/2,1-(1-conf.level)/2))
    return(list(median.1=median.1,median.2=median.2,median.diff=med.diff,hpd.diff=hpd.diff,quantile.diff=quantile.diff))
}

build.x.mat<-function(model,data.x.mat){
tt<-terms(model)
Terms<-delete.response(tt)
 m <- model.frame(Terms, data=data.x.mat, na.action = na.omit, 
            xlev = model$xlevels)
X <- model.matrix(Terms, m, contrasts.arg = model$contrasts)
return(X)
}


##################
###Regime Types
#########

#fmd1-Left Regime
data.x.mat.fml.left<-data.frame(femeale_exe1=0, prop_female=median(defense[rownames(model.matrix(model1.br)),"prop_female"]),
labor_full=median(defense[rownames(model.matrix(model1.br)),"labor_full"]),
mil_spend=median(defense[rownames(model.matrix(model1.br)),"mil_spend"]),
time=1:21,
left_exe1=1,
fmd1=1, 
deaths_any=0,
dictator1=0,  
involve=1,
wcombat=0,
oecd=0) 

yfmlleft<-cumprod(1-predict(model1.br, data.x.mat.fml.left,type="response"))
 yfmlleft[(c(1,5,10,15, 20))]

plot(cumprod(1-predict(model1.br, data.x.mat.fml.left,type="response")),type="l",lty=2,lwd=2,xlab="Time",ylab="Survival Probability",col="gray48") 

#fmd1-Right Regime
data.x.mat.fml.right<-data.frame(femeale_exe1=0, prop_female=median(defense[rownames(model.matrix(model1.br)),"prop_female"]),
labor_full=median(defense[rownames(model.matrix(model1.br)),"labor_full"]),
mil_spend=median(defense[rownames(model.matrix(model1.br)),"mil_spend"]),
time=1:21,
left_exe1=0,
fmd1=1, 
deaths_any=0,
dictator1=0,  
involve=1,
wcombat=0,
oecd=0) 


yfmlright <-cumprod(1-predict(model1.br, data.x.mat.fml.right,type="response"))
 yfmlright[(c(1,5,10,15, 20))]


points(cumprod(1-predict(model1.br, data.x.mat.fml.right,type="response")),type="l",lty=3,lwd=2, col="gray60") 


#Dictator
data.x.mat.dic<-data.frame(femeale_exe1=0, prop_female=median(defense[rownames(model.matrix(model1.br)),"prop_female"]),
labor_full=median(defense[rownames(model.matrix(model1.br)),"labor_full"]),
mil_spend=median(defense[rownames(model.matrix(model1.br)),"mil_spend"]),
time=1:21,
left_exe1=0,
fmd1=0, 
deaths_any=0,
dictator1=1,  
involve=1,
wcombat=0,
oecd=0) 

x.mat.model<-build.x.mat(model=model1.br,data.x.mat.dic) 

ydic<-cumprod(1-predict(model1.br,data.x.mat.dic,type="response"))
 ydic[(c(1,10,15, 20))]

points(cumprod(1-predict(model1.br, data.x.mat.dic,type="response")),type="l",lty=1,lwd=2, col=1) 

#legend(1, .6, c("Military Dictatorship", "Non-Left Former Mil. Dict.", "Left Former Mil. Dict."), col=c("black","gray60", "gray8"), lty=c(1,3,2),lwd=c(2,2,2))


#OECD

data.x.mat.oecd<-data.frame(femeale_exe1=0, prop_female=median(defense[rownames(model.matrix(model1.br)),"prop_female"]),
labor_full=median(defense[rownames(model.matrix(model1.br)),"labor_full"]),
mil_spend=median(defense[rownames(model.matrix(model1.br)),"mil_spend"]),
time=1:21,
left_exe1=0,
fmd1=0, 
deaths_any=0,
dictator1=0,  
involve=1,
wcombat=0,
oecd=1) 

x.mat.model<-build.x.mat(model= model1.br,data.x.mat.oecd) 

ydic<-cumprod(1-predict(model1.br,data.x.mat.oecd,type="response"))
points(cumprod(1-predict(model1.br, data.x.mat.oecd,type="response")),type="l",lty=4,lwd=2, col="black") 

legend(1, .6, c("Military Dictatorship", "Non-Left Former Mil. Dict.", "Adv. Indust. Dem.","Left Former Mil. Dict."), 
                col=c("1","gray60","1","gray48"), lty=c(1,3,4,2), lwd=c(2,2,2,2))
   

 ###########               
 ######TIME
 ###########

 quadfun<-function(t,m,a,b,c){
 	a*(t-m)^2+b*(t-m)+c
 }
a<-  -0.010141513 #Time Squared
b<-0.069155806 #Time
c<- -2.689350009 #Intercept
t<-c(1:21)

plot(quadfun(t,11,a,b,c)~t)
plot(exp(quadfun(t,11,a,b,c))~t)
plot(exp(quadfun(t,11,a,b,0))~t)


###################
###SUPPLEMENTARY INFO
#############

###################
###Supply of Prospective Female Cabinet Ministers
###################


#% Female Cabinet Ministers

summary(model.femin<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1 +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+wcombat+labor_full+self_appoint1+pwomenmin_int,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#High-Prestige Female Cabinet Minister

summary(model.hiprofile<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1+hiprofile,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))


#Suffrage

summary(model.suffrage <-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1+ suffrage,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Years Since Women Eligible to Serve in Political Office

summary(model.office_w <-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1+ office_w,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Years Since Women First Elected to Political Office

summary(model.elect_w <-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1+ elect_w,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

###################
###Recovery from Conflict
###################

#Post-Conflict Recovery

summary(model.rec<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+wcombat+labor_full+self_appoint1+ recovery,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Post-Conflict Recovery x Left

 summary(model.recleft<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
fmd1+
 as.numeric(involve)*
 log(mil_spend)+wcombat+labor_full+self_appoint1+  left_exe1*recovery,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Internal Post-Conflict Recovery

summary(model.recint<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+wcombat+labor_full+self_appoint1+ internal_recovery,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Internal Post-Conflict Recovery x Left

summary(model.recintleft<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
fmd1+
 as.numeric(involve)*
 log(mil_spend)+wcombat+labor_full+self_appoint1+  left_exe1*internal_recovery,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))


###################
###Cabinet Size
###################

summary(modelnum.min<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1+sizeofcabinet,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

###################
###Regimes Types
###################


#Democracy


summary(model.regime<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+wcombat+labor_full+self_appoint1+dictator1+dems,data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

# Former Civilian Dictatorships

summary(model.regime.fcd<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+wcombat+labor_full+self_appoint1+fcd1 ,data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#GDP (PPP)

summary(model.gdp<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+labor_full+wcombat+self_appoint1+dems+log(gdpppp),
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Presidential versus Parliamentary Systems

summary(model.pres1<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1+
 dictator1+dics+parliamentary_c1+mixed_c1,data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Excluding Botswana

summary(model.botswana<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1,
data=defense[defense$sname!="Botswana"&defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Coding Botswana as a Democracy

summary(model.botswana.dem <-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+ botswana.dem +
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Unified Governments

summary(model.unified2<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1+ unified,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Political Party Competitiveness

summary(model.comp<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1+herfgov,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))



###################
###Importance of Def. Ministry
###################

#CINC Score

summary(defense$cinc_ts[which(defense$AFFL_lead1==0&defense$nodef==0)])

summary(model.cinc<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
prop_female +femeale_exe1+
  left_exe1*fmd1+
  as.numeric(involve)+
  oecd+labor_full+wcombat+self_appoint1+log(cinc_ts),
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Conscription

summary(model.recruit1<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
prop_female +femeale_exe1+
  left_exe1*fmd1+
  as.numeric(involve)+
  oecd+labor_full+wcombat+recruit,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

#Military Officers

summary(model.mildef<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)+oecd+wcombat+labor_full+self_appoint1+mildef,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

###################
#Quotas
###################

summary(model.quota2<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female*quota +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))


###################
#Interstate Relations
###################

#EU

summary(model.eu<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full+wcombat+self_appoint1+eu,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],x=TRUE))

###################
#Conflict Measure
###################

summary(model.both.conflicts<-glm(FFL_lead1~I(time-11)+I((time-11)^2)+
 both.conflicts +dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+
 oecd+labor_full+wcombat+self_appoint1,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],family=binomial(link=logit),x=TRUE))

###################
#Self-Appointments
#############

summary(model2<-logistf(FFL_lead1~I(time-11)+I((time-11)^2)+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+
 oecd+labor_full+wcombat+self_appoint1,
data=defense[defense$AFFL_lead1==0&defense$nodef==0,],family=binomial(link=cloglog),x=TRUE))


###################
#Women's Appointment to Other High-Prestige Posts
#############

summary(model.hipi<-glm(FFHiPi~I(time-11)+
defense_appoint+
 prop_female +femeale_exe1+
left_exe1+
oecd,
data=defense[defense$AFFHiPi==0,],x=TRUE))

###################
#Women in Combat
#############

summary(model.consequence<-glm(wcombat~defense_serve+
 deaths_any+dictator1+
  prop_female +femeale_exe1+
 left_exe1*fmd1+
 as.numeric(involve)*
 log(mil_spend)+oecd+labor_full,
data=defense[defense$nodef==0,],family=binomial(link=logit)))

